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Abstract. We argue that the evolution of coloured partons into colour-singlet hadrons has approximate fac- 
torization into an extended parton-shower phase and a colour-singlet resonance-pole phase. The amplitude 
for the conversion of colour connected partons into hadrons necessarily resembles Regge-pole amplitudes 
since qq resonance amplitudes and Regge-pole amplitudes are related by duality. A 'Regge-cascade' fac- 
torization property of the A'^-point Veneziano amplitude provides further justification of this protocol. 
This latter factorization property, in turn, allows the construction of general multi-hadron amplitudes in 
amplitude-squared factorized form from (1 — » 2) link amplitudes. We suggest an algorithm with cascade- 
decay configuration, ordered in the transverse momentum, suitable for Monte-Carlo simulation. We make 
a simple implementation of this procedure in Herwig-|— |-, obtaining some improvement to the description 
of the event-shape distributions at LEP. 



PACS. 12.38.Aw General properties of QCD 
- 13.87.Fh Fragmentation into hadrons 
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1 Introduction 

By 'hadronization', we loosely mean the final phase in the 
process of jet formation, where coloured partons at the 
end of the parton shower turn into colour-singlet hadrons. 

The hadronization phase is presumably not completely 
separable from the perturbative phase, but approxima- 
tions can be made, and we can for instance set the per- 
turbative coupling to zero in the hadronization phase and 
vice versa. This is, in effect, the approach ordinarily taken 
in the Monte-Carlo event generators . 

A better approximation may be to allow the perturba- 
tive coupling to extend into the hadronization phase, and 
continue the perturbative evolution down to zero ^3^. 

The universality, i.e., that it does not depend on, for 
instance, where the other partons are, of the cut-off in 
the first case and the coupling in the second case can be 
justified in the Gribov confinement picture Although 
the complete treatment of the gluon Green's function is 
lacking at present, for the quark Green's function, Gri- 
bov's method indicates that its behaviour is governed by 
a universal equation containing both the gluonic semi- 
perturbative and the long-distance super-critical contri- 
butions. We can then separate out the semi-perturbative 
contribution by means of the effective coupling procedure, 
so that the remaining dynamics, of hadronization, would 
be a predominantly colour-singlet interaction, mediated 
by resonances and poles. In other words, the coloured par- 
tons that remain at the end of the extended parton-shower 
phase turn into hadrons by an interaction which is effec- 
tively colour-singlet. 



Colour preconfinement |SE] dictates that at the end of 
the parton shower, the colour-connected parton pairs, i.e., 
the colour dipoles, have a mass spectrum with a charac- 
teristic scale of a few times the parton-shower cut-off, and 
this mass scale normally comes out to be about 1 GeV. 
However, this is violated in low-p-r jets 0. 

As is the case in the Monte Carlo event generators, let 
us consider these colour-preconfined units as the starting 
point of hadronization. The colour-connected partons then 
exchange objects that are effectively colour-singlet to turn 
into hadrons. 

A justification for the resemblance of the confining dy- 
namics with colour-singlet exchange is in duality, i.e., the 
observation that the summation over resonance states re- 
produces the dynamical behaviour characteristic of Regge 
poles, and vice versa 




Fig. 1. A link in the cascade-decay chain. 



In this paper, we concretize this statement by observ- 
ing that the explicitly dual A^-point Veneziano (i.e., the 
open bosonic string) amplitude |10llll| satisfies an am- 
plitude squared factorization property that corresponds 
to cascade decay, where each vertex, shown in fig. ^ has 
Regge-like angular behaviour, except in the final decay. 
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Long-distance dynamics can then be treated as a se- 
ries of Regge-hke decay. This is Uterally true when the 
decaying unit is heavy, but because of semi-local duality 
|12| . Regge theory often remains a good approximation 
even in the low-energy region, or the low-mass region of 
the colour-preconfined units, where it is not formally sup- 
posed to be applicable, particularly after the resonances 
and dips have been integrated over. 

Combining this result with the above statement of ap- 
proximate factorization of semi-perturbative and hadroni- 
zation phases, we have a picture of hadronization where 
partons evolve via the parton shower and hadrons evolve 
via the Regge cascade. The exchange in the hadronization 
phase normally consists of the meson trajectories, but the 
baryons and even the pomeron trajectories are also al- 
lowed. 

The formahsm is developed between sees. |21 and 
There is ambiguity as to whether, and to what ex- 
tent, the forced g ^ qq splitting at the end of the parton 
shower, which is adopted by HERWIG ^ more or less for 
the sake of convenience, occurs in reality. In ref. we ar- 
gued that there may be physical origin to the low-energy 
enhancement of the g qq splitting, due to the string 
tension 'pulling apart' the colour octet. Another possi- 
bility would be that the gluon acquires a pole mass and 
hence can decay into two light quarks over a finite time. 
Even in this case, some of the gluons may remain intact, 
as all gluons do in PYTHIA 2 . Since finite-time effects 
(gluon decay) are not completely separable from infinite- 
time effects (colour-singlet interaction), the latter effect 
cannot be neglected. In sec. |H1 we discuss the hadronic 
observables that are potentially sensitive to this aspect of 
hadronization. 

In sec. |5| we examine the implementation of a simpi- 
fied procedure based on the preceding discussions in Her- 
wig-|--f . We compute a number of jet observables at the Z° 
pole and compare against the Hcrwig-t--f and LEP num- 
bers. 



2 The (2 2) Regge amplitude 

The amplitude for (2 2) scattering with Mandelstam 
variables s and t, omitting the signature factor, is: 



hadronization. in units of GeV, is then: 



A^mr{£-a{t)){~s/sor^'l 



(1) 



P{t) is the coupling factor. Either (3{t) or P{t)r {£ - a{t)) 
is often taken to be constant. a{t) is the ^-channel trajec- 
tory and £ is the spin of its lowest-lying member, sq is the 
Regge characteristic scale. We identify s with the dipole 
mass squared later on. 

The leading fiavour-singlet trajectory is the the near- 
degenerate p/ll>/ f /a family: 



a{t) ^ao + a't^ 0.5 + 0.9t, t=\. 



(2) 



t is measured in GeV^. When -^i becomes large, above 
10 GeV or so, we should also consider the pomeron con- 
tribution. The typical transverse momentum generated by 



{kT, 



hadronization 



l/v/0.91og(s/so). 



(3) 



For s « lOso, which can occur when, for example, — 1 
GeV and -y/s ~ 3 GeV, we have a transverse momentum of 
0.6 GeV. This gives a measure of the extent of the violation 
of local parton-hadron duality due to the colour-singlet 
phase. 

For lighter dipoles, the transverse momentum becomes 
greater. However, the transverse momentum cannot ex- 
ceed a half of the dipole mass. From eqn. we see that 
the turn-over should occur near s = sq. 

A natural choice of sq is obtained by comparing against 
the Veneziano model. Corresponding to eqn. (^), we have 
a Veneziano amplitude: 



A = (3 



r {£ ~ a{t)) r {£ - a{s)) 
r {£ ~ a{s) - a{t)) ■ 



(4) 



In the Regge limit, by applying the Stirling factorial ap- 
proximation, we recover eqn. Q with the extra constraints 
/3(t) = const, and sq = l/o;', where a' is the slope. For 
a' — 0.9 as in eqn. 10), we have so = 1.1 GeV^. Using 
the same a' in the s- and i-channels is justified since a' 
corresponds to the inverse mesonic string tension, which is 
physically, although not necessarily in Regge phenomenol- 
ogy, a universal constant. Sq = 1-1 GeV^ implies that the 
typical transverse momentum generated in hadronization 
is at most ^0.5 GeV, even for the lighter dipoles. 

For the sake of comparison, the HERWIG cluster mass 
cut-off has the default value of 3.5 GeV so that the typ- 
ical cluster mass is about a half of this, and so the typi- 
cal transverse energy is about 0.9 GeV. Clusters are the 
colour-connected units that form the seed of hadroniza- 
tion in HERWIG. PYTHIA has the kx distribution gen- 
erated artificially by a double Gaussian distribution, and 
the width a of the primary Gaussian distribution has the 
default value of 0.36 GeV. This does not imply that had- 
ronization in HERWIG is harder than in PYTHIA, as will 
be demonstrated by a simulation in sec.|Hl One of the rea- 
sons is that we have so far neglected the hadron masses. 

The ratio of the yield of hadron pairs that require heav- 
ier flavour exchange in the hadronization phase, to the 
yield of the states that only require the exchange of light 
flavours, is given as a function of their invariant mass Mj 
as: 

2 ^2(ao-0.5) 



2 

inv 



oc 



(5) 



Here ao is the intercept of the exchanged trajectory, and 
is less than 0.5. We have squared the amplitude to obtain 
the probability. We have assumed that the effect due to 
the difference in the slope a' of the two trajectories can 
be ignored. 

A possible practical application of the above formula 
would be in baryon pair production in jets. The largest 
ao for baryons is 0.0, corresponding to one of the A/N 
trajectories ^Sl- The relative baryonic yield is therefore 
proportional to the inverse of the invariant mass squared. 
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The physics of three-body baryonic decay of iJ-mesons 
[T^ is subject to the same consideration. We can under- 
stand the enhancement of the three-body decays to two- 
body decays as the suppression of the baryon-exchange 
when the baryon pair mass is large. This view is supported 
by the so-cahed 'threshold effect', i.e., the tendency that 
the baryon-antibaryon pair is formed with small invariant 
mass. A similar phenomenon is seen in the production of 
pp in low-virtuality 77 collision at Belle |15| . where ad- 
ditional mesons often accompany pp, particularly when 
sufficiently above the threshold [T^. 



3 Factorization of the A^-point amplitude 




N-1 



N-2 



Fig. 2. The A^-point amplitude. 

The iV-point Veneziano amplitude corresponds to the 
scattering of open bosonic string shown in fig. [21 This is 
often expressed in a cyclically symmetric form but 
for our purpose, the formulation of ref. |17j is more con- 
venient. We have: 



N-2 1 



i=2 \ k=j 



(6) 



The approximations made in the above equation, in terms 
of the trajectory, is a{t) — —^^ + t where jj, is the mass 
of the external legs. This approximation will be removed 
in the factorized formula to be derived later on. g^~'^ is 
the coupling factor, ti are the square of the momenta that 
flow in between the i'th and the i + I'th legs, defined by: 



N 



(7) 



The multi-Regge limit is obtained in eqn. © by the 
approximation that Ui are small. However, when consid- 
ering the application to hadronization, this multi-Regge 
formula is not very convenient, for two reasons. The first 
is that resonances are not incorporated. The second is that 
it is difficult to construct an iterative evolution algorithm 



based on this equation, that starts from two incoming ob- 
jects. 

Let us derive a different limit of eqn. ((HJ, corresponding 
to the cascade-decay configuration illustrated in fig. 13 




N-2 
N-1 

Fig. 3. Cascade-decay picture of the A'^-point amplitude. 

In eqn. @, we consider a resonance in the Vih. and 
Z + I'th legs. We write: 

2p( • pi+i + /i^ = {pi + Pi+iY - = "(si) =n- e. (8) 

n is an integer which corresponds to the maximum spin 
in this channel. The integral is dominated by the region 
ui ^ I. If « 1, we find that: 

lim A(l,...,l-l,Ll + l,...,N) 

^ A{1,... ,1-1,1 + 1,. ..,N) 

r-l 



X g lim^y dui (I-u;)-i-"+^ 



ccay ■ 



(9) 



The amplitude thus factorizes into the production part 
and the decay part. This is distinct from the usual res- 
onance factorization property of the iV-point amplitude 
prij . since the latter is in general not amplitude-squared 
factorizable and is therefore of limited use. 

The decay part is just the Euler beta function. With 
the usual analytization convention, this becomes: 



hm A, 
£^±0 



r{-n + e)ri~a{ti)) 

decay -5 _ ^ ^ ^) • U^j 



The form of the amplitude is identical with the 4-point 
Veneziano amplitude excepting the difference in the cou- 
pling factor. If n is large, the angular distribution has 
Regge behaviour. Resonance leads to the factorization in 
eqn. but the Regge angular distribution is indepen- 
dent of resonance. In general, so long as the mass of the 
decaying system is large, one can consider the dynamics 
to be dominated either by resonances, with u/ w 1, or by 
poles, with ui ~ 0. 

We can therefore consider a 'Regge cascade' iterative 
procedure, in which we have a series of two-body decays. 
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each of which except the final ones exhibiting Regge angu- 
lar distribution. We can avoid intermediate configurations 
that are manifestly non-Regge by requiring that for each 
decay, s is large and the masses of the two decay products 
are small. This can be achieved by the choice of a sen- 
sible ordering variable, related to the mass. This will be 
discussed in sec. 

Let us consider the production amplitude in eqn. (jSJ. 
This is as given by eqn. ©, but the mass term needs care. 

We consider a link in the cascade-decay chain that has 
resonant daughters. Denoting the decay by [0 12] as 
shown in fig. Q the decay distribution is given by: 



^docay(0 ^12) 



— q lim 



(11) 



This is as before. However, eqn. ||SJ) needs to be modified 
to: 



2pi ■P2+ 



/ I \2 2 2,2 



(12) 



In terms of the 4-point Veneziano amplitude, we have the 
replacement: 



a{s) a{s) - a{si) - a{s2), 
so that the decay amplitude is modified to: 
r{-a{s))r{-a{t)) 



(13) 



decay 
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r{-a{s) - a{t)) 

^ rj^ajs) + ajsi) + a{s2))r{-a{t)) 

^ r{-a{s) + a{si) + a{s2) - a{t)) 

The Regge limit of the amplitude is given by: 



(14) 



^decay(0 ^ 12) ~ gr{-a{t)){-a{s) + a{si) + a{s2)r^'l 

(15) 

For simplicity, we have omitted other permutations of 
external particles. The introduction of some of these so- 
called 'twisted' terms result in the signature factors: 



(_«(,))"(*) ^ («(,))"(*) 



lie 



-i7ra=^(t) 



(16) 



The sign between the two terms, i.e., the signature, de- 
pends on the nature of the trajectory being exchanged. 



4 The density function 

We now turn to developing a practical algorithm based on 
the factorization of eqn. and the decay amplitude of 
eqn. H15|) . We first calculate the mass distribution for the 
daughters 1 and 2 in fig. H by the optical theorem. 

We temporarily introduce the decay width F in order 
to keep track of the phase space factors. The two-body 



decay width of an object with mass ^/s, expressed in terms 
of t, is given by: 



dr{0 12) 1 



rl^(0^12) 



dt 167rs3/2 
The density function is defined by generalizing this to 
dr{0 12) 1 



(17) 



dtdsids2 IGtts'^/^ 



\A{Q^l2)\'p{si)p{s2). (18) 



We evaluate p{s) by the optical theorem. Starting from 
the total decay width, which is in general given by: 



r = 



1 1 



^-^ImKO^X^O)], 



(19) 



we obtain: 



p(si) = ^Ini[A(l->X^l) 



(20) 



The amplitude A(\ X ^ 1) can be estimated as: 
r(-a(si) - a(0)) 



(21) 

so that after substituting — a' (3, we have: 

Pisi) « ^r(-a(0))(a(5i))"(°) sin(^a(0)). (22) 

Let us absorb the F function and the phase factor in /?. 
There are several possibilities for estimating this coupling 
coefficient. For instance, an order-of-magnitude estimation 
can be obtained by comparing against the tt'^tt" total 
cross section. We have: 

(T*°*(12 -^X) = ilm^(12 -> 12)|t=o 

w ^(a's)"(o) = a'/3(a's)"(o)-\ (23) 
s 

From ref. 9 , using coupling factorization, we have: 



(T 



tot, Roggoon/ 4- 



{tt+tt-) 



a{TT'^p)a{'K p) 

aipp) 
27.56 X 36.02 

56.08 * 



-0.4525 



Hence: 
so that: 



/3 w 20mbGeV'^ w 50, 
50a' 



(mb)(24) 



(25) 



27r 



^0.5475 ^ 7„0.55 



(26) 



p is measured in units of GeV"^ and s is measured in 
units of GeV^. p, in principle, includes the non-continuum, 
resonance, contribution as well as the other trajectories. 

The density thus obtained does not agree with the one 
obtained by integrating over two-body decays. The con- 
tribution to the density function from two-body decay is: 

P(s)o^l2 = / ^^^j^^^^ |^dccay|^p(si)p(s2). (27) 
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This does not agree with eqn. H26|l . This situation is fa- 
mihar from perturbation theory. The single emission of a 
gluon from a parton gives rise to an infrared-divergent con- 
tribution to the cross-section. This is resolved by adding 
together the virtual corrections of the same order. The 
Sudakov form factor gives an all-order expression for the 
probability of emission (or no emission), and this effec- 
tively sums, up to an infrared cutoff, the divergent parts 
of the real and virtual diagrams together. 

Making the correspondence with the perturbative case, 
eqn. (|20|l is the all-order sum. Eqn. (|27|l corresponds to 
the single-emission cross section. The 'splitting function' 
is the ratio of the two: 



f splitting (<Z')d<^ 



1 dtdsids2 
p{s) IBtt^s 



|^docayrp(si)p(s2X-28) 



d'P is the phase space element as before. Writing down 
the Sudakov form factor requires the choice of a suitable 
evolution variable. 



5 The ordering variable 

The factors to be considered in proposing the ordering 
variable for cascade decay are: 

1. The algorithm generates all configurations and avoids 
double counting. 

2. The ordering should be physical, i.e., there should be 
Regge behaviour at every vertex except the last decay. 
Vertexes where the daughter masses are comparable 
with the mass of the mother should be avoided where 
possible. 

3. It should be local, i.e., the decay in one branch of the 
tree should not depend on the decay of the other. 

In view of the above, it would seem that ti defined in 
eqn. © may be suitable, as any collection of ti has unique 
ordering so long as there are no degenerate subsets. Thus 
the first point is satisfied. We may further improve this or- 
dering and say that the ordering is in each branch of the 
tree, so that in order to generate the decay somewhere in 
one branch, it is not necessary to look up the values of ti in 
other parts of the tree. This then satisfies the locality con- 
dition. If we require that \ti\ are ordered in the decreasing 
order, the decay with the largest daughter masses would 
occur first, so that the requirement of Regge behaviour at 
every vertex would be approximately satisfied. 

However, the form of eqn. (|28|l suggests that this may 
not be the best choice insofar as the ease of event gen- 
eration is concerned. A better choice would involve the 
squared masses si and S2- We therefore add the condi- 
tions: 

4. The ordering variable is expressible in terms of the 
masses of the daughters. 

5. An approximate solution is acceptable if the approxi- 
mation has a physical ground. 



This suggests the quantity: 

SlS2 



\t\. 



(29) 



Unlike ti, Si depend on the configuration of the tree, and 
so does, to some extent, the above quantity. This would 
sometimes lead to the violation of the uniqueness of or- 
dering. 

The quantity defined above is similar in form to the 
transverse momentum, so that we also choose to call this 
quantity k'^. This is proportional to the transverse mo- 
mentum of the emission that would be required for this 
cascade decay, had the decay been caused by the emission 
of a qq pair. We also define rapidity y as: 



1 



y = 7:log(s2/si)- 



From the form of p{s), we have: 

p{si)p{s2) _ 
Pis) ~ 



p(4). 



(30) 



(31) 



Using this relation, eqn. (|28|l simplifies to: 

^'spiitti„g(<Z')d^ = -^dk^dy l^dccayl^ p{kT). (32) 

Since |t|min ~ fcr ^'^'^ ^decay, from eqn. H15|l . only has 
mild dependence on si and S2, we conclude that the dis- 
tribution of the splitting function is almost flat in y. 

The Sudakov form factor, that is, the probability of no 
decay in between two specified phase space boundaries, is 
given in general by: 



Z\(f?) = exp 



(33) 



Using the simplified expression of eqn. ()32|l and further 
imposing the approximations: 



\t\r 



Adocay = 9r{-a{t)){-a{s)r^'\ (34) 



we obtain the estimation: 
Z\(fc|n, s) fa exp 

21og(a(s)) 



dQ^dy 



167r2 



(«(.)) 



2a{- 



^PiQ' 



(35) 



This formula then gives the approximate cascade-decay 
evolution. From this equation, we see that the dominant 
dynamical factor is in the exponent (q;(s))^". The inte- 
gration over y approximately yields logs, which cancels 
against the corresponding expression in the denominator. 
Absorbing the other dependence in /3 and after inte- 
gration, we obtain: 



Z\(fc^, s) fa exp 



(/?/100) 
log a's 



(36) 
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This can be contrasted with the perturbative expres- 
sion [H]: 



A{y) = exp 



-Tyrdy — - — 



(37) 



Comparing eqns. (|35|l and (|37|l . we may obtain an effective 
and non- universal 'non-perturbative as'- 



6 Algorithm for cascade decay 

The proposed algorithm for generating multi-hadron final 
states is, in outline: 

1. Terminate the parton shower, for instance, by means 
of a universal couphng. In the course of this process, 
some or all of the gluons become qq. 

2. Form colour-singlet units from colour-connected par- 
tons. Depending on whether they involve gluons or not, 
these would respectively correspond to kinked strings 
|2] and clusters Kinked strings may, for instance, 
be treated dipole-by-dipole. One first chooses a colour- 
connected dipole in the kinked string and, during or 
after its hadronization, turn to the resolution of the re- 
maining colour by re-interacting with the other colour- 
connected parton(s). 

3. The dipole hadronization proceeds by cascade decay. 
The evolution variable is fcy, defined by eqn. 1)291) . An 
initial value is set for this, at about s/4 where s is the 
mass squared of the dipole. The simplified form of the 
Sudakov form factor is given by eqn. 1)36(1 . 

4. Rapidity is generated either as a flat distribution, or 
by using, for instance, eqn. H32|) . The masses of the two 
daughters are calculated from fc^ and y. t is generated 

according to l^dccayl^- 

5. For each daughter, the initial value for the evolution 
variable is the of the decay that just took place. The 
t-channel 4-momentum is recorded so that the kine- 
matics for the next decay can be evaluated without 
having to sum over the momenta of other branches. 

6. At some stage, the daughters are identified with physi- 
cal resonances. This can be incorporated into the den- 
sity function p. 

7. The process is repeated until every branch either has 
been identified with a physical resonance or its has 
reached zero. The latter would correspond to a decay 
into two ground-state mesons. 

Although our starting point involved an explicit for- 
mulation of the multi-particle amplitude, the above pro- 
cedure is more general and depends only on the principles 
underlying the amplitude, namely duality, factorisability 
and Regge behaviour. The Regge amplitudes are gener- 
alizable to non-scalar ground-state hadrons, and it is a 
simple matter to include flavour. 



7 A Regge fragmentation function 

As a special case of the cascade-decay tree configuration, 
let us consider the 'fragmentation-function' configuration, 





N-1 N-2 



Fig. 4. Fragmentation-function picture of the iV-point ampli- 
tude. 



shown in fig. 01 This potentially involves cases with very 
heavy daughters so that, for instance, eqn. ()15|l may fail, 
but we proceed by treating it as an approximation. In 
place of eqn. ()18|l . we have: 



dr{0 12) 
dtds2 



1 



167rs3/2 



|A(0^12)|Vs2). (38) 



The first daughter is the physical, ground state, hadron 
that splits off. We neglect the mass of this hadron. The 
splitting function is: 



„ , ^ , p{s2) dsn dt , , ,2 

Psplitting(f )d<? = ^ V 16^ ^"^"-'^ ■ 



(39) 



We denote the energy fraction of the branching particle, 
i.e., particle 1, by z, so that S2 = s(l — z). Our approxi- 
mation is valid when z w 1. We obtain: 

Pspiitting(<^)d^ = (1 - z)0-55dz|^(a(sz))2"(*). (40) 

The evolution variable this time is z, so we integrate over 
t. For a(0) = 0.55, we obtain: 

iDTT'^ 21og(a(sz)) 
The approximate Regge fragmentation function is: 

fiz) = P{z) « z^-^{l - z)°-55 X 0(1). (42) 



8 Signature of gluon splitting in the 
three-body final state 

As discussed in the Introduction, the amount of g ^ qq 
forced splitting is an ambiguity in our picture of hadroni- 
zation. In this section, we consider the hadronic observ- 
ables in low-energy low-multiplicity events, three-body in 
particular, that are sensitive to this aspect of hadroniza- 
tion. Since HERWIG and PYTHIA represent two extreme 
parametrizations of this splitting, we carry out a simula- 
tion using these generators. Both generators are expected 
to perform badly with the default parameter set, as few- 
body final states are not what the generators are designed 
for. 

We consider the study of exclusive three-body pri- 
mary hadron production in charm events at BELLE^. 

^ We thank A. Chen for suggesting the study of these events. 
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In both generators, we adopt the default parameter sets. 
The centre-of-mass energy is 10.58 GeV corresponding to 
BELLE, and we select the charm pair production events, 
with the matrix element correction turned on so that the 
hard e"'"e~ —> ccg configuration is generated according 
to the perturbative matrix element. We generate 10000 
events in each case. 

We turn ofi^ the decay of hadrons in order that the 
generated hadrons are 'primary'. The generated final state 
corresponds to the reconstructed few-body events with no 
further resonances among the final state particles. 



Final state 


HERWIG 


PYTHIA 


1 hadron 








2 hadrons 


313 


269 


3 hadrons 


359 


1293 


4 hadrons 


2810 


2795 


5 hadrons 


955 


3340 


6 hadrons 


3828 


1699 


7 hadrons 


483 


494 


8 hadrons 


1150 


96 


> 8 hadrons 


102 


14 


Even no. hadrons 
Odd no. hadrons 
All 


8162 
1838 
10000 


4860 
5140 
10000 



Table 1. The number of primary hadrons, generated using 
HERWIG and PYTHIA. 



We first present a table of the number of the pri- 
mary hadrons in tab. |H| HERWIG, based primarily on 
the isotropic two-body decay algorithm, normally leads 
to an even number of hadrons in the final state. The ex- 
ception to the rule occurs when one of the clusters is too 
light to decay into two hadrons. In this case, the cluster is 
identified with the lightest meson with the corresponding 
quantum numbers. This predominance of even-numbered 
multiplicity is considered to be an artifact I !l , but since 
this phenomenon is due to the forced g qq splitting, 
it is possible that an effect of this nature may occur in 
reality to some extent. 

Both generators predict a typical multiplicity of about 
four to six primary hadrons. 

We now look at the kinematic distribution of the three- 
hadron final state, rejecting the rare events containing 
baryons. When doing so, we propose to make use of a 
two-dimensional scatter plot in rapidity y and the loga- 
rithm log(fcT) of the transverse momentum, both being 
measured against the 'jet' direction. We choose this set 
of observables in order that perturbative emissions have 
an almost flat distribution. In the soft limit, we have the 
DGLAP gluon emission probability 

a-Pcmission = "7^"^ ■ (43) 

C is the colour factor which, for gluon emission from a 
quark line, is Cf- We often adopt sa k'^. 



We carry out the simulation at the BELLE centre-of- 
mass energy of 10.58 GeV. For both the charm quark mass 
and the charmed meson mass, we choose 2 GeV. 

In the three-body case, it is more useful to replace 
the rapidity and Ut by quantities that are more directly 
measurable. For massless particles in the soft limit, for 
the reaction [0 123], we write in terms of the invariant 
masses: 

= " , 2/ = log {M23/M13) . (44) 

-'"123 

Particle 3 is soft, and can be identified with the non-charm 
meson. We adopt these definitions even for massive final 
states. One advantage of doing so is that a plot in the 
(log kx, y) plane is the Dalitz plot in the (Mfg, MI3) plane 
on the log scale, rotated by 45 degrees. The resonances, 
which form horizontal and vertical lines on the Dalitz plot, 
form straight lines with slope ±1 on this plot. 

log(kj/Ecm) vs. y 

-0-5 P ' Pert. ' ^ ^ 
.PYTHIA 

-1.5 - ■ Y 

-!■ * JK' - , i'^:hk 

'+*+r^ * * * ^ > fl. 

-2.5 - 

-3 - 



Fig. 5. The log(fcT/\/s) ^s. y scatter plot for charm events. 
We show the results of simulation based on the perturbative 
three-body matrix element with modified as (red), PYTHIA 
(green) and HERWIG (blue). 

In fig.El we plot the results of HERWIG and PYTHIA, 
along with the results of simulation based on a pertur- 
bative three-body matrix element for e^e~ —> ccg using 
the modified as of ref. [3- The choice of the renormal- 
ization scale for the perturbative result is = k'^. The 
gluon is identified with the soft particle 3. The peak of 
as is at 0.575 GeV. The number of generated events is 
1000 for the perturbative case and 10000 for HERWIG 
and PYTHIA. However, because of the low probability of 
three-body events, the number of accepted events is much 
less for both generators. 

On the vertical axis, a transverse momentum of 1 GeV 
corresponds to w —2.3, and 0.5 GeV corresponds to sa —3. 
The perturbative distribution is almost flat in rapidity as 
discussed above, with some modification due to the finite 
mass. 
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Turning to the behaviour of the Monte Carlo event 
generators, we see that the behaviour of the two generators 
is quite distinct. 

The PYTHIA distribution is easier to understand. The 
band at near 1 GeV, corresponding to —2 on the vertical 
axis, is understood as a combination of the perturbative 
gluon emission and the artificially generated px, related to 
the Gaussian width parameters PAR J (21) - PARJ(24). 

The HERWIG distribution has a complicated struc- 
ture. The structure resembling faint lines at ±45 degrees 
indicate the typical mass of the cluster that decays into 
one charmed and one non-charmed meson. The distribu- 
tion is in general softer than the PYTHIA distribution 
because the extra particle comes from the decay of a clus- 
ter and most of the energy in this decay is taken up by 
the charmed meson. There is little remnant of pertur- 
bative gluon emission since all of the gluons have split 
into qq and these recombine with the colour-connected 
quark/antiquark to form clusters. 

As mentioned at the beginning of this section, the 
Monte-Carlo predictions should not be trusted for few- 
body exclusive final states. Nevertheless, the experimen- 
tal investigation of such configurations can shed light not 
only on the areas in which the generators could be im- 
proved but also on the mechanism of hadronization as a 
whole. On the other hand, there would be considerable 
and, seemingly ^Hl, not insurmountable difficulty asso- 
ciated with the experimental reconstruction of primary 
hadrons within the plethora of decay products. 



9 Implementation in Herwig++ 

We now turn to the implementation of the algorithm de- 
veloped in the previous sections. We choose to adopt Her- 
wig-|--|- [2n| as the platform. The hadronization algorithm 
of Herwig-h-h is based on that of HERWIG. 
Our implementation entails: 

1. Modified daughter mass distribution in cluster decay, 
generated according to eqn. H3t)|) . which can be in- 
verted analytically. We impose ordering in kx, and the 
minimum daughter mass is the sum of the constituent 
quark masses with the corresponding flavour. Clusters 
that are too light to be decayed in the default Her- 
wig-| — h implementation are not decayed. 

2. Option of 'cluster rotation', or the Regge smearing of 
the directions of the 3-momenta of the daughter clus- 
ters, as well as the partons making up the daughter 
clusters, with respect to the partons making up the 
parent cluster. This has the distribution cx (a's)^" *. 
Here t is the Mandelstam variable for the (2 2) sub- 
process involving, in the initial state, parent cluster 
partons, and in the final state, the daughter clusters. 

3. Option of pomeron inclusion. By imposing the approx- 
imation 2q!pqjjjj,j.qj^ — aRcggcom resulting Sudakov 
form factor is still analytically invertible. 

4. Option of Regge flavour and baryon generation taking 
place after cluster splitting, with weight cx (a's)^"^"^ 
and adjustable constant multiplicative factor for each 
flavour and diquark. 



In eqn. H36(l . we use the parametrization: 

(3h = 13/100^1, a' = 0.9GeV-^ a(0) = 0.5. (45) 

Out of the three parameters, a' and a{0) are fixed by 
Regge phenomenology. The remaining parameter, (3h, is 
not well determined, but the physical results typically de- 
pend only on the logarithm of fin- This is because the 
typical kx of cluster decay is the dominant factor in the 
decay kinematics. The Herwig-|--|- parameter CLPOW be- 
comes redundant. We start by keeping the values of the 
other parameters the same as in the default set. 

All our simulation results are for e'^e~ interaction at 
the Z° pole, corresponding to the LEP centre-of-mass en- 
ergy. We generate 100,000 events in each case. Our simu- 
lation results are compared against the default set of ex- 
perimental data in Herwig-|--|-. Both soft and hard 
matrix-element corrections are switched on. 

We first show our simulation result for 1 — T in fig. El 
After including the cluster-rotation effect, the agreement 
with data is comparable with that of Herwig-f-l-, and is 
better in the region of low 1 — T, corresponding to pencil- 
like events. We found similar improvement in most other 
event-shape observables, with the exception of oblateness. 
For spherical events, not shown in the figure, the results 
are similar to those obtained using Herwig-|--|-. 



1-T 



10 




With cluster rotation 
No cluster rotation 
Hw++ default 
DELPHI 96 



0.02 



0.04 



0.06 



0.08 



0.1 



Fig. 6. The distribution of 1 — T, with and without the cluster 
rotation effect. 



In addition to the event shapes, we examined the be- 
haviour of the single particle and the identified particle 
distributions. For the former, the results tend to be better 
than the Herwig-|--|- numbers for rapidity, but worse for 
transverse momenta with respect to thrust and sphericity 
axes. As an example of the latter, fig. [7| shows the in- 
plane transverse momentum with respect to thrust axis. 
The transverse momentum distributions improve when we 
omit the cluster rotation, although doing so would ruin the 
agreement with the event shapes. One speculative possibil- 
ity would be in the over-generation of transverse momenta 
in our combined Regge-HerwigH — h implementation. 
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In-plane transverse momentum w.r.t. thrust axis 



100 



10 



0.1 



0.01 



0.001 



1e-04 



With cluster rotation 
No cluster rotation 
Hw++ default 
DELPHI 96 



Fig. 




7. 



10 



12 



14 



The distribution of in-plane transverse momentum 
with respect to thrust axis, with and without the cluster rota- 
tion effect. 



For the identified particle distributions, the results are 
similar to those of Herwig-|--|-. The distribution of scaled 
momentum of charged particles, shown in fig. |S1 shows 
deficit compared with the experimental data both in the 
very large- and small- regions where Xp stands for the 
scaled momentum. For the large- region, we find better 
agreement when the sample is restricted to fight-, charm- 
or bottom-quarks, so that this is possibly not a serious 
deficit. In the small- Xp region, the agreement is improved 



1000 



100 



With cluster rotation 
No cluster rotation 
Hw++ default 
SLD 03 




0.01 



Fig. 8. The distribution of the scaled momentum of charged 
particles, with and without the cluster rotation effect. 



by lowering the HerwigH — h shower cut-off parameter S. We 
found that the flavour dependance of the shower cut-off is 
milder in our approach. HcrwigH — h uses the parametriza- 
tion: 

Qg^iS- 0.3m,)/2.3, (46) 

where Qg is a gluon virtuality cut-off and S is constant, 
but this yields too much radiation from 6-quarks and too 



little from light quarks in our approach. The choice of 
S = 2.3 for 6-quarks and 6 = 1.7 for c-quarks was found 
to be more appropriate. An alternative possibility is to 
introduce the pomeron. We found that the substitution: 



(as) 



2qe 



0.02(a^ 



(47) 

reproduces the observed charged-particle multiplicity. Do- 
ing so, however, leads to the deterioration in the descrip- 
tion of the very large-Xp region. On the other hand, this 
deterioration in the large- Xp region makes the agreement 
with the thrust-like event-shape distribution better in the 
pencil-like region. 

The distribution of the proton was found to be in poor 
agreement with the data and this is similar to the case 
of Herwig-f-t-, but the agreement is better in the case of 
proton production in charm events. This is shown in fig. El 
This seems to imply that we may have an incomplete de- 
scription of either the shower cut-off or hadronization in 
our Regge-Herwig++ amalgamate approach, in the case 
of the light parton from the hard processs. 



100 



The scaled momentum of protons for charm events 



With cluster rotation 
Hw++ default 
SLD 03 



The scaled momentum of charged particles for all events 




0.01 



0.001 



Fig. 9. The distribution of the scaled momentum of protons 
for charm events. 



For the last item of our implementation, namely that 
of flavour generation during cluster decay, we found no 
visible effect on any of the distributions. 



10 Conclusions 

We argued that the dynamical conversion of partons to 
hadrons can be effectively factorized into two phases. The 
first is an extended perturbative phase based on an univer- 
sal infrared coupling. The second is a long-distance hadro- 
nization phase mediated by colour-singlet resonance-pole 
dynamics. 

From the consideration of the Veneziano A^-point am- 
plitude, we argued that, essentially because of duality, 
amplitude-squared factorization is applicable to the lat- 
ter phase as well as the former phase. This factorization 
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has the form of cascade decay, where each decay except 
the last one in the chain has Regge angular behaviour. 

Since both phases are factorizable, we can describe the 
overall fragmentation process by merging together the two 
factorized contributions. 

By generalizing the factorization of the A/^-point am- 
plitude, we proposed a framework for Monte-Carlo event 
generation. We derived an approximate Regge fragmenta- 
tion function as a special case. 

Our results are applicable generally to multiple hadron 
production from pole-resonance dynamics. 

We made a simple implementation of our hadroniza- 
tion algorithm in HerwigH — h. We found that the event- 
shapes tend to be better compared with the original pro- 
cedure. On the other hand, the results are mixed in the 
case of single-particle and identified-particle distributions. 
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